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0\ \ ABSTRACT 

nJ , This is a review of recent work on the dynamic response of Josephson junction arrays 

(_P ■ driven by dc and ac currents. The arrays are modeled by the resistively shunted Joseph- 

[^^ ' son junction model, appropriate for proximity effect junctions, including self-induced 

^_^ ' magnetic fields as well as disorder. The relevance of the self-induced fields is measured 

xj , as a function of a parameter K ^ Xj^/a, with A^ the London penetration depth of the 

Oa ■ arrays, and a the lattice spacing. The transition from Type II (k > 1) to Type I (k < 1) 

behavior is studied in detail. We compare the results for models with self, self+nearest- 

C^ ' neighbor, and full inductance matrices. In the K ^ OO limit, we find that when the 

S , initial state has at least one vortex-antivortex pair, after a characteristic transient time 

I 1 these vortices unbind and radiate other vortices. These radiated vortices settle into 

'"^ ' a parity- broken, time-periodic, axisymmetric Coherent vortex state (ACVS), char- 

ed ' acterized by alternate rows of positive and negative vortices lying along a tilted axis. 

Q , The ACVS produces subharmonic steps in the current voltage (IV) characteristics, typ- 

^ ■ ical of giant Shapiro steps. For finite K we find that the IV's show subharmonic giant 

Shapiro steps, even at zero external magnetic field. We find that these subharmonic 

steps are produced by a whole family of coherent vortex oscillating patterns, with their 

t^^ , structure changing as a function of K. In general, we find that these patterns are due 

'%_] ■ to a break down of translational invariance produced, for example, by disorder or anti- 

j^ ' symmetric edge-fields. The zero field case results are in good qualitative agreement with 

experiments in Nb-Au-Nb arrays. 

1. INTRODUCTION 

There has been considerable contemporary interest in the study of the dynamic 
response of two-dimensional Josephson junction arrays (JJA) subjected to external 
probes 1-27 -pj^jg attention has been due in part to recent advances in photolitho- 
graphic fabrication of these arrays with tailor-made properties. The dynamics of the 
JJA are described in terms of a large set of coupled, driven, non-linear differential 
equations. When the current drive has dc-|-ac components, novel non-equilibrium 
stationary coherent vortex states may appear in the arrays. Experimentally, giant 
Shapiro steps (GSS) have been observed in the IV characteristics of proximity-effect 
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JJA in zcrocl and rational magnetic field fi'ustrations cl. The fi'ustration, / = <I>/$Oj 
is defined as the average applied magnetic fiux per plaquettc, $, measured in units 
of the fiux quantum $0 — h/2e. Coherent oscillations of ground state field-induced 
vortices are believed to be responsible for the existence of the fractional GSS when 
/ = p/?) with p and q relative primes Q. This interpretation was successfully verified 
in numerical sipjulation modelling of the arrays bv_a resistively shunted junction 
(RSJ) model Icl, as well as from analytic studies El. There have also been some 
experimental E£l and theoretical studies O'HEJ of the relevance of the geometry of 
the arrays and the direction of the input current on the generation of fractional 
GSS. Moreover, the experiments in GSS have encouraged investigations of JJA as 
coherent radiation sources EJ. 

Later experiments on Nb-Au-Nb EjOO and Nb-Cu-Nbi''') JJA also showed 
half-integer GSS (and also some evidence for higher order subharmonic steps) in 
zero magnetic field. These steps can not be explained within the context of the 
non-inductive RSJ model, which m-pyed successful in providing an understanding 
of the / — p/q Shapiro steps. Q'Q'LlO'tj In this review we will show that these extra 
subharmonic steps can be produced by breaking the translational invariance in the 
arrays. 

For example, we found that the additiDXLoLdisorder can in fact induce extra 
half-integer steps in zero magnetic field ll2't2lE3 within the RSJ model, even in 
the limit when no self-field effects are included. These extra steps were found 
to be related to a coherent oscillation of current-nucleated vortices. Previously, 
the nucleation of vortex-antivortex pairs by the defects, induced by an applied 
dc current, was discussed by Xia and Leath Ell. In our studies we added an ac 
component to the dc current. The addition of the ac current changes the vortex 
dynamics in a fundamental way: It leads to the formation of a novel, far from 
equilibrium, axisymmetric coherent vortex state (ACVS). The ACVS corresponds 
to an oscillating pattern of tilted positive and negative vortex rows which produce 
the extra half-integer GSS (see Figs 6 and 7). The relevance of the disorder is to act 
as nucleating centers of vortex pairs, thus providing the relevant initial conditions 
in the dynamics to generate the ACVS. However, as we shall discuss in this review, 
what matters in generating the ACVS is to have a mechanism to produce vortex- 
antivortex pairs (VAP) in the initial conditions since then the dynamics are such 
that the ACVS is generated in most cases. We will discuss in more detail the 
properties of the ACVS in Section 3. 

Most of the experiments in GSS have been carried nut in proYir nity effect arrays, 
that have strong temperature-dependent critical currents IjOEJO. In this case self- 
induced magnetic field (SIMF) effects, which also break translational invariance due 
to the current induced edge-fields, can h£,of relevance for the interpretation of the 
experimental data as was found in Ref. 113. We have recently developed a dynamic 
model to study the SIMF effects in Josephson junction arraysE3EjO. When applied 
to the study of dc-|-ac driven JJA, we found that the fractional GSS observed at 
nonzero magnetic fields are affected by SIMF effects. There are two extreme regimes 
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of interest as a function of k = A^/a in which the underlying microscopic coherent 
vortex states responsible for the fractional GSS are qualitatively different. Here A^ 
is the London penetration depth and a the lattice spacing. 

Also, in the case of zero external field we found subharmonic GSS, produced by 
a farnjly of oscillating coherent vortex states, with different structures as a function 
of K t3. This connects with the ACVS which turns out to be the k ^ oo limit of 
this family of vortex states. 

In this review we present a brief recap of our main results pertaining to the 
generation of coherent vortex states, mostly in zero external magnetic field. In this 
way we present a unified view of our understanding of the physics of subharmonic 
Giant Shapiro Steps at / = 0, with zero and finite screening current effects. In our 
initial calculationscj we considered the cases where the inductance matrix had diag- 
onal and nearest neighbor components. This approximation has been signi, 
improved by including the full inductance matrix in the analysis of Refs I 
Following the general idea of this improvement we have redone some of our pre- 
vious calculations including the full inductance matrix and we include a critical 
comparison between the results of the different approximations. 

2. DYNAMICS OF NON-INDUCTIVE JJA 

2.1. The Josephson effect 

We start by briefly reviewing the essentials of the Josephson effect. A Josephson 
junction is made of two small superconductors separated by a thin film of non- 
superconducting material. The two basic equations that describe the physics of the 
junction are 

/j = /osin(02-ei), (1) 

that gives the supercurrent that flows between the two superconductors, and the 
voltage drop between them given by 

V.||(«.-«0. (2) 

Here 9i , 02 are the phases of the Ginzburg-Landau order parameter in supercon- 
ductors 1 and 2, respectively. The critical current /q is the maximum current that 
can flow through the junction, and $o = h/2e is the superconducting quantum 
of flux. In the presence of a magnetic field the phase difference is replaced by its 
gauge invariant form, {62 — ^1) ^ (^2 — 6*1 — -^ /^ A- dl), where \/ x A = B is the 
magnetic field. In real junctions there is always some dissipation. This dissipative 
effect has been successfully modeled, in the case of proximity effect junctions, by 
the McCumber-Steward resistively shunted junction (RSJ) model c3. In the RSJ 
model, for a dc current biased junction, the total current / that flows in parallel 
with the ideal Josephson current is 

/ = /,/ + Yin = /o sin( A0) -f ^ ^^ , (3) 
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where we have used Eqs. (nl) and (0) and TZ is the RSJ shunt resistance. When / is 
time-independent, Eq.(H) can be solved analytically. In this case, the time averaged 
voltage {V) as a function of /, which defines the IV curve of the junction, is given 
by {V) = when / < /q and {V) = Tly/ 1"^ - I^ when / > /o. 

When the junction is driven by a time-periodic current / = Idc + ^ac sin(27ri^i), 
the {V) vs Idc curve shows plateaus at the quantized voltages 

(K)=«^, n = 1,2,3,... . (4) 

These are thii Shapiro steps that have allowed very precise measurements of the 
voltage unit cj. The central question addressed in this paper is what happens when 
we couple a large number of Josephson junctions to an array that is then driven by 
dc-f-ac currents. 

2.2. Josephson junction arrays 

A Josephson junction array ( J JA) is made of a.n NxN network of superconduct- 
ing islands connected by Josephson currents Id'H. For example, the square arrays of 
Benz et al. Q were made of 1000 x 1000 Nb-Cu-Nb proximity effect junctions with 
a lattice constant of a = 10 fim. A schematic representation of a JJA driven by an 
external current is shown in Fig. 1. 

To model the dynamical behavior of this system we extend the RSJ model to 
a square JJA network 0. The current /^(r), along the fx direction between the 
superconducting islands at sites r and r + /i, with /i = ix,ey, is given by 

I^{r) = lo sm{A,e{r) - A,{r)) + ^ | {A,0{r) - A,{r)) + r,,{r, t) , (5) 

where A^6{r) — 9{r+^)~9{r). The external magnetic field produces the frustration 
f = ■^^- , that measures the average number of fiux quanta per unit cell, and defined 

by 

27r/ = Ax{r) + Ay{r + e^) - A^ir + iy) - Ay{r) = A,, x A^,{r), (6) 

where the link variable A^(r) — ^ JJ '^ A ■ dl. Here, for the moment, we are 
neglecting the screening currents by assuming that Afj^{r) is fully determined by the 
external magnetic field H. This assumption is correct whenever Al ^ Na. For the 
same reason, -^ can be dropped from Eq. (0) since the magnetic field is constant 
in time. In Sec. 4 we shall present the full analysis including screening current 
effects. The effects of temperature are included by adding the Gaussian random 
variable 77^ (r, t) to the equations of motion with covariance 

{ri,{r)v,'{r'))^^ryS,y6{t-t'). (7) 

Eq. (ph together with Kirchoff's current conservation law, 

A^ • I^ir) = I,{r) - Ix{r - x) + Iy{r) - Iy{r - y) = I^^^r), (8) 
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valid at each node, fully define the evolution of the phase 6{r, t) as a function of 
time. Here I'^^^{r) denotes the external current injected at site r. The explicit 
expression for ^ derived from these equations is u 

^^ = -^EG(''''^')U^^*(^,i)-A,,-[/osin(A^%',t)-A^(r')) + ^p(r',i)]}, 

r' 

(9) 
with G{r, r') being the two-dimensional lattice Green function which depends on the 
boundary conditions chosen. Eq.(ph defines the set of coupled nonlinear dynamical 
equations studied in this paper. In Benz's experiment U this would represent a set 
of 10^ non-linear coupled oscillators driven by external currents. 

Of course, solving this large set of coupled non-linear equations analytically is 
out of the question. We use then an efhcient algorithm to study Eq.(H) numerically. 
In our simulations the external currents are injected uniformly along the y-direction 
with /^^*(y = 0) = I{t) and zero elsewhere. The potentials at y = Ny are fixed 
to zero by setting the phases 0{y = Ny -I- 1) = (this boundary condition removes 
a singularity present in the Green function matrix). We choose periodic boundary 
conditions (PBC) along the x-direction in most calculations, but we have also used 
free end boundary conditions (FBC) for comparison purposes. 

A direct numerical evaluation of Eq. (^) grows as iV'*, which is too slow and 
limits the sizes of the lattices that can be simulated. Overcoming this restriction 
turns out to be essential for the formation of the ACVS state. We have used in- 
stead a very efficient method ll3't3E3 based on the use of the Fast Fourier Transform 
(FFT). Specifically, this involves doing a FFT of Eq. (J^) along the x-direction and 
then solving the resulting tridiagonal matrix equation along the y-direction. This 
algorithm grows as N"^ log N , which is appreciably faster than the direct method. 
Eikmans and van Himbergen Q were the first ones to use the FFT along both the 
X and y directions to study JJA. Our method shows an improvement over theirs 
of about 30%. The time integration is carried out using a fixed step fourth order 
Runge-Kutta (RK) method. Furthermore, at finite temperatures we use an exten- 
sion of the second order RK method for stochastic differential equations developed 
by Helfand and GreensideEd. Typical integration steps used were Atz^o = 0.01 — 0.1, 
with characteristic Josephson frequency vq = 21310. _ 

2.3. Periodic square arrays and giant Shapiro steps 

For ordered arrays, Eq.(|9|) is invariant against the transformations 

f^f + n and / ^ -/. (10) 

with n an integer. In this case it is enough to analyze the properties of the model 
in the interval / = [0, 1/2]. The phase vortex excitations that can appear in the 
JJA are defined by 

g [A^0{r)~A^ir)] = 2n{n{R)-f), (11) 

R 
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where the sum is over the plaquette R and the gauge invariant phase differences, 
^A'^('') — ^/i(^), are restricted to the interval [— 7r,7r]. Therefore, n{R) is an integer 
that gives the vorticity at plaquette R. In zero field and at finite temperatures 
vortices plajj a central role in triggering the Berezinskii-Kosterlitz-Thouless phase 
transition eH. For fractional frustrations / = p/q, with p and q relative primes, 
the ground states of the airays consist of superlattices of field-induced vortices 
with unit cell of size g x g l2I. These states lead to non-monotonic behavior of the 
magnetoresistance and Tc{f) studied both experimentally and theoreticallytlo'El. 

We can calculate the experimentally measurable IV characteristics from the 
solutions to Eq. (^). The calculations give the time-averaged total voltage drop {V) 
per row in the array as a function of the applied dc current. For periodic arrays, 
the / = dynamic response of the model is equivalent to the superposition of 
N individual Josephson junctions along the direction of the current^ for the total 
current flows along the y-direction. The model behavior is that of N^ rows with Ny 
junctions in series. Thus, when the lac = 0, the IV characteristics are simply given 
by the one-junction result multiplied by Ny. The same is true in a JJA when / = 
and I{t) = Idc + lac sin(27r:^ai), leading to giant Shapiro steps at voltages u: 



nhv 
"27 



yn = Ny^^, n= 1,2,3,... . (12) 



The Ny factor is important in the possible practical applications of JJA. The sym- 
metry that allows us to separate the iV^; x Ny array into Nx independent rows in 
series is broken when the external magnetic field is nonzero. No exact analytic 
solution to Eq. (g) that gixes the IV characteristics is known in this case. It was 
first found experimentally eI that there are fractional GSS at voltages 

K.g = ^^, n,g= 1,2,3,... , (13) 

q 2e 

when / — p/q. These fractional steps were explained as due to the collej 
oscillations between different / = p/q vortex ground state configurations Erl 
The vortex-superlattice oscillates in synchrony with the ac current with frequency 
h'/q. We have in fact seen these ground state oscillations in animations produced 
with solutions to Eq. (|^). 

In the zero field case one could expect that by breaking translational invariance, 
which allows the reduction of the 2-D array into effectively 1-D rows, one may get 
extra subharmonic GSS even in zero magnetic field. One can break translational 
invariance by simply adding disorder to the array. This is the topic of the next 
section. Later, in Sec. 5, we will see that current induced magnetic fields also break 
translational invariance. 

2.4. Defect nucleated vortex pairs. 

We introduce disorder in a square JJA by displacing the lattice sites radially 
away from their periodic positions by a distance S < 1, and then choosing the 
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angle from a uniform probability distribution in [0, 2tt]. We assume that the critical 
currents and shunt resistances are all the same and equal to Iq and 7?.'s, respectively. 
This means that the disorder is only evident when the link variables A^(r) ^ 0. 
Adding one defect to the array affects only the four v4^ 's connecting the bonds linked 
to the displaced lattice site. Varying / can be seen as increasing the deformation or 
disorder in the lattice, since the A^'s are proportional to the external field. When 
we take f — n an integer, because of the local periodicity of the equations of motion 
with respect to /, only the plaquettes around the defect will be modified by the 
presence of the field. The n value chosen will represent the "strength" of the defect. 
As we shall discuss later the specific nature of the disorder will not be important 
in producing an ACVS. For example, the ACVS also appears when we cut bonds 
with / = 0. 

Leath and Xia Ell studied the response of a JJA with a rectangular defect and 
in the presence of a constant dc current. They found that there is a tendency to 
create a vortex- antivortex pair (VAP) pinned to the defect for Idc 7^ 0. When the 
external current is larger than the critical current of the array, /c, the Lorentz force 
acting on the VAP is enough to break them apart and they can move away from 
the defect. As they move along the direction perpendicular to the external current, 
they produce a dissipative Faraday voltage that leads to a nonzero resistance in the 
array. A larger number of VAP's is generated when Idc ^ Ic, leading to rather 
complex vortex motions in the array. 

In Fig. 2 we show an example of this situation for the kind of topological defect 
we defined above. We show the vortex distribution n{R) for a 40 x 40 array, which 
was the typical size studied, although we will discuss results for larger lattices as 
well. For I < Ic we see a VAP pinned at the defect. When I > Ic, vortices and 
antivorticcs move away from the defect as found in Ref. t3. 

3. AXISYMMETRIC COHERENT VORTEX STATES (ACVS) 

3.1. IV characteristics 

In this section we discuss the IV characteristics of JJA with the type of defect 
introduced in the previous section and driven by a current 

/""* = /dc + /acSin(2^z/t). (14) 

Let us start with a JJA with only one defect. Later we will discuss what happens 
when the JJA has more defects. 

In Fig. 3(a) we show the IV characteristics for a periodic JJA. In calculating 
the IV characteristics as we vary Idc, we take the final phase configurations of the 
preceding current value as initial conditions for the next one. In Fig. 3(a) we see 
the expected integer GSS (Eq. [ij), which are identical to those obtained with one 
Josephson junction, except for a factor of Ny. In Fig. 3(b) we show the results for 
the IV of a 40 X 40 square lattice with one defect located at the center of the array. 
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We can clearly distinguish the new half-integer steps in the IV curve at voltages 

i,/2 



n N,,hiy 
K^/2=2^; " = 1,2,... , (15) 



in addition to the expected rounding of the usual integer steps due to the presence 
of disorder. It turns out that the precise location of the defect on the lattice has 
no effect on the nature of the final stationary state nor on the appearance of the 
half-integer steps. The half-integer steps are found to be hysteretic. This is seen 
in Fig. 3 (c) where we show a blow-up of the 1/2-step region, but for an 80 x 80 
lattice. For the 1/2-step we can define two critical currents, /+ when increasing 
the current, and /" when decreasing it. As usual, the size of the hysteresis loop 
is history-dependent. For / > /+, the 1/2-step is reversible, and it can be reached 
independent of the initial conditions. For I^ < I < I^ the response depends 
directly on the initial conditions. Note the appearance of a 2/3-step at 

which was not clearly evident in the 40 x 40 lattice calculation. The 2/3-step is also 
found to be hysteretic. We will discuss this step in more detail in Sec. 3.3. 

An important property of the half-integer steps is that they are not clearly 
visible for lattice sizes smaller than abouLTS x 18. We have then carried out a 
finite size analysis of the 1/2-step width E3£3 and found that there is a minimum 
critical size Nc ~ 16, above which we start to see indications of the existence of the 
1/2-step. The 1/2-step width reaches an asymptotic plateau for lattices larger than 
about 32 X 32. We believe that the need to have a minimum lattice size to see the 
1/2-step is related to the specific nature of the ACVS structures that needed to be 
stable in the non-equilibrium stationary state. We discuss this point further in the 
next section. 

We have also calculated the TV's for JJA with free boundary conditions. We 
find that the subharmonic steps are also present, but the lattice sizes needed to see 
the ACVS are still larger. This is easy to understand since the boundary has more 
of a perturbing effect in smaller lattices. As the lattice size increases the shape of 

^ubharmonic steps approaches the sharpness of the ones found when using PBC 
Therefore, we conclude that the boundary conditions are not essential for the 
existence of the ACVS. 

Once we had identified the subharmonic steps due to the presence of defects we 
studied their quantitative properties. Following the approach used for one junction, 
we began the characterization of the half-integer Shapiro steps by studying their 
width-dependence on the frequency i>, the ac current lac, the amount of disorder 
and temperature. In Fig. 4(a) we show the dependence of the 1/2-step width as a 
function of i/ with lac = ^o. Here we note that there is a frequency window for which 
the 1/2-step is clearly visible and has a well defined maximum at fmax ~ 0.22 i/q li . 
The width of the frequency window is Ai/ = j/+ — i/_ ^^ 0.0431^0 — 0.4i/o. Outside 
this frequency window the step width is essentially zero. 
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In the one junction case the step width as a function of lac goes hke A/„ = 
2/o'/n(^^f^)j with J„ the Bessel function of integer order Eil. In Fig. 4(b) we show 
the oscillatory behavior of A/1/2 ^ a function of lac, for i^ = 0.1 vq, and with one 
defect. The oscillations do not appear to fit the expected Bessel function form but 
the data is not good enough to make a definitive statement. Even more important is 
the fact that we appear to find a critical value for lac above which the step appears 

We need to consider temperature effects to make comparisons with the exper- 
imental data. We have then studied the stability of the 1/2-step as a function of 
temperature E3£3c3. We found that the 1/2-step gets rounded as temperature in- 
creases and the hysteresis step width is reduced. The 1/2-step tends to disappear 
for a relatively low temperature T > 0.04$o^o/27rfcB, which is small compared to 
the BKT critical temperature (~ 0(1)) Eil. This is consistent with the fact that the 
width of the 1/2-step is small when compared to the critical current of the array. 

We also found that when we added more disorder, by increase the number of 
lattice sites displaced from their periodic positions in the array, the half-steps are 
still there but their hysteretic properties get quantitatively modifiedll3£j. Even 
when all the lattice sites are deformed, the 1/2-integer step is still there, but it is 
more rounded and with the hysteresis loop reduced in size. We conclude that a 
completely distorted lattice is not enough to destroy the extra 1/2-steps. Howeveii 
a study of the half-integer steps as a function of the strength of the disorder Ej 
shows that there is a critical value above which the 1/2-steps do disappear. This 
value coincides with the disorder strength above which there is a plastic flow of 
vortices in dc driven JJA t3. 

3.2. Coherent Vortex oscillating states 

In the previous section we described the conditions under which extra subhar- 
monic steps appear in the IV characteristics in the k = 00 limit. Here we want 
to discuss the physical origin of these steps based on a microscopic analysis of the 
vortex dynamics. 

A useful quantity to look at is the spectral function 

S{v)^ lim I- r V{t)e'^'"''dt\^ (17) 

as a function of frequency ly. We calculated S{p) and found, as would ha|ffi been 
expected, that it shows resonances at frequencies 7t,z//2 for the 1/2-steps tBiEi. The 
resonances at the 1/2-step indicate that there is an underlying coherent vortex 
oscillatory state with twice the period of the ac current. These resonances only 
exist for currents right on the steps, fractional and integer, while for other values 
of Idc there is a broad band spectrum in S{i'). 

As mentioned in Section 2.3 the presence of fractional giant Shapiro steps, found 
when / = p/q, is explained in terms of coherent vortex oscillations between different 
lattice ground states. At zero field and at T = there are no vortices in a periodic 
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array. However, as we showed in Sec. 2.4 for the dc case, the presence of defects can 
produce current-induced vortices. Therefore, to further study the vortex dynamics 
we calculated the total absolute number of vortices in the array 

A^aW = ^|n(i?,t)|, (18) 

R 

as well as the total vorticity 

NT{t)=^7i{R,t)^N+ + N_, (19) 



where n{R,t) was defined in Eq. (|ll|). We started by choosing Idc within the 
reversible part of the 1/2-step. Typical results for Na{t) and Nxit) in the single- 
defect case are shown in Fig. 5. From this figure we can identify two characteristic 
times td and tAcvs, the dissociation time and the ACVS time, respectively. For 
t < td there is a vortex pair (as in the dc case. Fig. 2) but with its ± polarity 
oscillating in phase with the current. The size of the VAP increases with time due 
to the oscillating Lorcntz force felt by the vortices. At td the vortex-dipole breaks up 
and the vortices move away from the defect. After that time the number of vortices 
increase steeply. There is a novel mechanism that generates vortex pairs from the 
moving ones. This kind of vortex "radiation" phenomenon appears to be catalyzed 
by the presence of the oscillating field due to the presence of the dc-t-ac current. 
The vortex motions in this intermediate phase do not follow a regular pattern. In 
the one-defect case the transient state has a left-right antisymmetry about an axis 
centered about the defect and along the current direction. This antisymmetry is 
reflected in Nt, which is exactly zero in this time interval. After the tAcvs time 
the antisymmetry is broken, and an axisymmetric coherent vortex state (ACVS) 
starts to form. This is seen in the Na(t) curve as periodic variation of the number 
of vortices, and in NT{t) as fluctuations about zero. 

In Fig. 6(a) we show the distribution of vortices and antivortices in the ACVS 
state. The vortex configurations are formed by rows of vortices of alternating sign 
lying along a tilted axis with an angle of about 27° measured clockwise with respect 
to the X-axis. The rows are separated along the x-axis by a distance of about Nx/2. 
In total there are Ny/2 vortices and A^j,/2 antivortices. This state oscillates in sign 
in phase with the ac current. After a time l/v the vortex rows interchange their 
signs, and after 2/^ they are back to the previous vortex state. Therefore, the 
vortex configurations change periodically in time with a period twice that of the ac 
current. During most of the time the vortex patterns correspond to one or the other 
possible ACVS configuration. The change of sign of the vortex rows happens in a 
very short time, corresponding to the peaks in the Na{t) curve of Fig. 5(a). At that 
instant in time the vortices and antivortices appear to cross from one row to the 
other. The collisions between them produce additional vortices and antivortices, 
which annihilate each other rapidly before a new ACVS of opposite sign is formed. 

The ACVS remains stable for the largest times considered (10"*/;/). To study 
the stability of the ACVS, we turned off the ac or the dc currents at times toff 
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as shown in Fig. 5. In both cases, after turning off one of the currents with the 
other one still on, the ACVS collapses into the initial VAP pinned state. The 
total vortex number decreases via direct vortex-antivortex annihilations. However, 
the process of annihilation is much slower when the ac current is still on while 
the dc current is turned off. This means that the ac current strongly affects the 
collisions between vortices by actually delaying the direct annihilation of vortices 
by antivortices. Within a 1/2-step the ACVS is also stable against small changes 
in the frequency and magnitude of the ac current and the frequency. 

A remarkable aspect of the ACVS state is that its structure, characterized by 
its angle, the distance between vortex rows and the number of vortices, appears 
to be a very robust non-equilibrium fixed point attractor of the dynamics. For 
example, we also have considered the cases when two, three, or all the lattice sites 
are disordered. We find that in these cases, although the quantitative values of the 
transients and the number of disorder-induced vortices is different, the stationary 
oscillatory state is essentially an ACVS. Moreover, there are parameter windows for 
the ac current, frequency, defect "strength" / and temperature, where the ACVS 
shows the same vortex pattern. The periodically oscillating vortex pattern is also 
the same for higher order n/2-steps with period (n + l)/^. However, we needed 
larger lattices, at least 64 x 64, to "fit in" an ACVS when the boundary conditions 
were changed to FBC, although the ACVS patterns appeared slightly "distorted" 
from the ones obtained using PEC. We conclude that boundaries can affect the 
coherence of the state, but for large enough lattices the ACVS forms and is stable. 

We note that the ACVS can form in the two possible broken-symmetry states 
with angular orientations either 27° or its complimentary 180° — 27°. Furthermore, 
we found that for large rectangular lattices Jjdien Ny >> N^) the two broken- 
symmetry states can coexist with each other 1130. 

At this point we wanted to know what essential ingredients are needed to nu- 
cleate an ACVS. A naive approach would suggest that the steady state is reached 
when the rate of vortex generation by the defect is equal to the rate of vortex an- 
nihilation. However, we found that most vortices were generated away from the 
defect during the transient time {td < t < tAcvs)- If disorder was essential in the 
formation of the ACVS, then when we switch off the disorder the state should col- 
lapse. However, the ACVS remained stable when we removed the defects by setting 
/ = 0. We also tried switching off the defects for times td < t < tAcvs and the 
ACVS is still stable. On the other hand, for times t < td the state collapses after 
we set / = 0. This suggested that once we have a large enough VAP, so that the 
vortices can not annihilate each other, the VAPs are then capable of generating a 
whole ACVS by themselves. This proves that what is essential in the formation of 
an ACVS is the "radiation" mechanism that produces and annihilates VAPs during 
their collisions. We also checked that if we start with a VAP in an ordered lattice, 
the ACVS is generated for the same Idc^ lac and v as when a defect was present. It 
appears then that: 
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Any physical mechanism (not only from disorder) that is capable of producing 
a vortex- antivortex pair sufficiently far apart, can nucleate an ACVS in a two- 
dimensional J J A. 

3.3. Higher order ACVS 

We already mentioned that wlien increasing the size of the lattice there are new 
subharnionic steps in the IV characteristics. This is shown in Fig. 3 (c), where an 
additional voltage step appears at 

for a lattice of size of 80 x 80. We also found that this step appears clearly after a 
minimum lattice size, which is about Nc ~ 70. We found that the time evolution of 
the vortex patterns responsible for the 2/3-step also has the general structure of an 
ACVS, but now oscillating with period 3/i/. This state becomes stable only after 
a very complex transient with a longer tAcvs than the one found in the 1/2-step 
case. In Fig. 6(b) we show the stationary vortex distributions for a current within 
the 2/3-step. It consists of rows of vortices (and antivortices) lying along a tilted 
axis with the same angle of about 27° with respect to the x-axis and separated 
by a distance of N^/'i along the x-direction. The total number of vortices is now 
|iVy. To have period 3 there is one row of vortices and two of antivortices. They 
interchange positions in synchrony with the ac-current. The total vorticity is still 
zero on the average, for the rows with vortices have double the density of those with 
antivortices. This implies that there may still be more possibilities for complicated 
ACVS-like oscillating patterns with the broken symmetry states in larger lattices. 
We note, however, that the step-width of the 2/3-step is smaller than that for the 
1/2-steps. One can speculate that the higher subharmonic states possibly present 
in larger lattices will have step- widths that will tend to zero as the size of the lattice 
tends to infinity. The IV characteristics could then have a Devil's staircase type of 
structure. 

3.4. Experimental observations of subharmonic GSS at / = 

In the previous sections we have seen that fractional giant Shapiro steps, pro- 
duced by oscillating collective vortex states, can be generated either by having a 
fractional frustration / ~ p/q in periodic arrays or by specific initial conditions that 
contain a vortex-antivortex pair (VAP) in the / = case. There we showed that the 
physical origin of the two types of collective vortex states is physically very different. 
These two scenarios for producing GFSS do not appear to be able to explain the 
lents that show Shapiro steps at fractions n/2 and n/3 in zero magnetic field 
3''^^ilLce these subharmonic steps were not observed in single Nb-Au-Nb 
junctionfOHEj, they must be due to the non-linear coupling between junctions in 
the two dimensional arrays. Therefore, they cannot be understood within the RSJ 
model given in Eq. (g) , because in the / = case the dynamics of a periodic array 
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is reduced to that of an effective single Josephson junction. However, as we dis- 
cussed in the previous section, the ACVS produces subharmonic steps in zero field 
within the RSJ model, provided there is at least one VAP in the initial conditions. 
Therefore, an ACVS could be a good candidate to explain the experimental results 
observed in Ref. EaliaO'El'^''^ . However, these experiments were done on (nominally) 
ordered arrays. More importantly, as it was recognized in the experimental papers, 
the critical currents in the SNS samples were large and therefore self-induced mag- 
netic fields, not included in the RSJ model, could play a role in explaining these 
extra steps [M3. 

The experiments in GSS have been carried out, for the most part, in oroximity 
effect arrays that have strong temperature-dependent critical currents la'Ej'E30. A 
carefully controlled narrow temperature range was studied experimentally in order 
to minimize the effect of the self-induced magnetic fields (SIMF) u . For example, in 
Benz's experiments the value of the London penetration depth quoted in reference 
4(b) is \l{Tc = 3.5K) = 280 (this estimate is obtained from using Pearl's formula, 
which is strictly valid in the continuum limit). Note that Xl/o, < N = 1000, 
which is in principle out of the regime of validity of the model given in Eq.(P) (e.g. 
Xl ^ Na). Xl decreased further when the temperature was lowered below Tc- 
Moreover, when going down to T ~ 2.5K, a zero-field 1/2-step was observed by 
Benz U'''\ with a value for A^ = 2.5a, which is certainly far from the region_^of 
applicability of the RSJ model. Also, in the experiments by H. C. Lee at al. lJ a 
value of Ai = 2.2a was quoted. 

3.5. Antisymmetric edge magnetic fields 

The discussion given above suggests that the SIMF can be essential for the un- 
derstanding of the experimental data. In fact, it was found in the experiments of 
RefJlj that the currents flowing in the array (at / = 0) produce an induced anti- 
symmetric magnetic field distribution in the arrays. The magnetic field is stronger 
at the edges of the array along the direction parallel to the external current, with 
magnitude Bedge ~ +MP^*/a'^ at one edge, and Bedge ~ —MP^^/a^ at the other. 
In the center of the array the magnitude of B tends to be small. Here M is the 
mutual inductance between plaquettes in the array. Note that these current in- 
duced magnetic fields, due to their antisymmetric nature, also break translational 
invariance in the J J A. 

A phenomenological way of trying to mimic the edge-fields, is to add them as 
boundary conditions to the non-inductive RSJ model. We have done numerical 
studies adding the edge-fields withflux $ — 7/'^^*$o at one edge of the array 
and $ = — 7/'^^*<I>o at the otherE3o. Again we found extra 1/2-steps in the IV 
characteristics for all 7 parameter values studied (0.05 < 7 < 2), just as we found 
before in arrays with defects. When looking at the animation of the vortex patterns 
responsible for these steps we found the same type of ACVS as discussed in the 
previous section. An specific example is shown in shown in Fig. 7. This pattern 
oscillates with period 2/;/ when the current is in the 1/2-step. Looking at the 
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transient that produces this ACVS, we found that the edge-fields produce VAPs, 
and then the dc+ac currents proceed to help nucleate the ACVS state, although 
the transient is longer and more complicated than in the case with defects. 

This result is in agreement with the conclusion we reached at the end of Sec. 3.2 
in that what is needed to produce an ACVS is simply a mechanism to initially 
nucleate VAPs in the array. These are the initial conditions that appear to lead to 
the ACVS. Note that in the calculations described here we used FBC. 

4. DYNAMICS OF INDUCTIVE JJA 

We just discussed how an ACVS can be produced by artificially imposing the 
edge magnetic fields at the boundaries of an array described by the RSJ equations of 
motion. However, the real antisymmetric magnetic fields measured experimentally 
must be produced dynamically from Faraday's law. In order to properly take into 
account screening current effects we must solve self-consistently the fiux and phase 
dynamical equations. This is the goal of this section. 

Here and in Sec. 5, we will describe our results from including the SIMF at 
different levels of approximation, together with a qualitative comparison to the 
experimentalresults. As mentioned in the introduction, our initial approximate 
calculationsE^ to the inductance matrix w£f6 recently superseded by including the 
full inductance matrix in the calculationsc3. Here we present our way of including 
the full inductance matrix into the calculation, which differs in the details from the 
approach followed by Phillips et al., but that it leads to similar physical results. We 
also present a comparison of our results obtained within the different inductance 
matrix approximations. 

4.1. Self-consistent JJA dynamical equations 

We start by noticing that the current along the bonds of the lattice, /^(r, t), is 
still given by Eq.(5). The difference comes only from the fact that in this case 

^0, 



dA^r,t) 



dt 

where the vector potential A^{r, t) is now related to the total magnetic flux $(i?, t) 
at plaquette R: 

^""^i^' ^' = A^ (r, t) ^Ay{r + e^,t)~ A^ {r + ey,t)- Ay (r, t) = A^ x A^{r, t) . (21) 

In this case we see that the total ^(i?, t) now depends on the flux generated by the 
external field $2; = Ha^ = f^o, plus the magnetic fiux induced by all the currents 
flowing in the array. Specifically we can write 

$(i?, t) = $,(i?) + J2 r(^, r, /x')/m' K> t) > (22) 

where T{R, r' , /i') is a matrix that explicitly depends on the geometries of the array 
and the junctions. It is convenient to write Eq. (E2|) only in terms of dual or 
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plaquette variables. From the current conservation condition given in Eq. (@), we 
have that 

I^ir,t) = J{R,t)-J{R-ey,t) 

Iyir,t) = J{R-e,„t)-J{R,t) + I,,t{t), 

with J{R,t) the plaquette's current, defined to be positive for currents flowing in 
the anticlockwise direction. We rewrite these equations in a short-hand notation, 

I^{r) = A^ X J{R) + 5^^yl,,t : (23) 

indicating that the external current is only applied along the y-direction. In terms 
of the plaquette variables Eq. ( p2|) can be written as, 

$(i?, t) = $,(i?) + ^ LiR, R')J{R', t) + E{R)h,t{t) (24) 

R' 

where 

L(i?,i?') = A^, xr(i?,r',/i') (25) 

is the inductance matrix of the array and we defined 

E{R)h,t{t) = ^r(i?,/,ey)/e,t(t) . (26) 

r' 

This term gives the magnetic flux induced by the applied external currents which 
is antisymmetric along the cc-direction. These magnetic fields have maximum am- 
plitude at the edges of the array and decrease towards its center, and thus they are 
called "edge magnetic fields" . For example, for a current sheet they are given by, 

„/„ „ N lina , .NtG — Rt, 
EiR,,Ry)^t^H "^^ " ), 

with R,j: e [a,(Nx — l)a]- This result is a good approximation for the actual value of 
E{R) in JJA l3. Note that the magnitude of E depends directly on the size of the 
array, and it decreases at the center of the lattice as lattice size grows. In the limit 
\^ - Rx\-^ N^a, E{Rx) « |f^(^ - i?:r) at the center of the array We note 
that there are other ways of separating and interpreting the different contributions 
in Eq.dl), as in Ref. H 

We note that to write down the set of dynamical equations for A^ (r, t) and 
6{r,t) we need to fix a gauge. In the Coulomb gauge, A^ • Ap(r, i) ~ 0, and from 
the combination of Eqs. (|), (|2^), (H) and (H) we obtain again Eq. (|) for 0ir,t), 
but now complemented by the following dynamical equation for the fiux 

d^(R,t) . ,. , ^ , . ^^ 

^ '' ^ nioA^xsin{A^9{r,t)-A^{r,t)) 



dt 



7^A^ 



J2 L-\R, R!) (*(i?', t) - *x - E{R')Iext) 



-e(i?,0(27) 
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with the noise function, 

C(i1■,i)=7^A^x?7^(r,i). 

The complete dynamics of the array is now governed by Eq. (0) and Eq. (|27 
These equations have to be solved for 6{r,t) and ^{R,t) self-consistently, where 
each one of these variables has its own characteristic frequency. We can define 
typical values for these frequencies from linearizing the Josephson term in Eq. (p7ffi3 
giving, vg = 2ttTZIc/^o = i^o, and v^ = TZ/^oa. To carry out an efficient numerical 
integration of the equations in time it is essential to take into account the relative 
magnitudes of the two relaxation times. In the extreme Type II regime i^^ ^ vg and 
the fast decaying variables are the fluxes while the 6''s are slow, with the opposite 
happening when ;/$ <C I'g- This situation is typical of 'stiff' problems in ordinary 
differential equations, which are notoriously difficult to treat analytically and even 
numerically, for they lead to singular perturbations c3. On the other hand, since 
the equations of motion are gauge invariant we can use this symmetry to find the 
most appropriate gauge to solve the problem. It turns out, however, that a fixed 
gauge does not allow us to efficiently solve the equations for all values of v,^/vg. In 
the v<j,/ve S> 1 limit, the extreme Type II regime, a convenient gauge to choose is 
the Coulomb gauge. We have implemented an algorithm that works in this case. 
Our discussion here will concentrate on the intermediate regime 0.05 < Vi^jvg < 20, 
where the stiffness problem i slesfLSf vere Furthermore, this is the regime considered 
in the experiments of Ref. OllaO, and it corresponds to the temperature range, 
[1.5 iC, 2.6 -ftT] in Benz's experiments eI'''''. In this v^/vg range it is more convenient 
to use the temporal gauge to rfficiently solve the equations. This gauge has been 
used extensively in the past cJ and more recently^ within the context of a JJA 
model of ceramic high temperature superconductors. This gauge entails replacing 

In terms of ^^ the Langevin dynamical equations of motion read, 
^^d^vW ^ _^A vyL-i(i?,i?')(A,x*Jr') + 27r/ + ^£;(i?')/e.*) 

R' ^ 

-lo sin 'J'^(r) + S^^ylext - Vt.{r, t) . (28) 

This system of equations describes the same physical dynamical evolution as Eqs. (Q) 
and (p7|), (note that in both gauges there are 2N'^ dynamical variables). Of course, 
we have to take free boundary conditions in order to allow the internal and external 
magnetic fields to relax to their correct stationary values. This implies that rela- 
tively large lattices have to be simulated to get the physically correct asymptotic 
behavior. 

4.2. Inductance matrix models 

A general analytic expression for the inductance matrix L{R, R'), valid for arbi- 
trary [R, R') and array geometry, is not known. However, we can learn a lot about 
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the main qualitative properties of the array from the asymptotics of L{R, R') and its 
general symmetries. We start by writing the standard definition of the T{R,r' , fi') 
matrix defined in Eq (pa) 



T{R,r\t,')^^——i ^^" :: ' T"7^ ' . (29) 

47r bRbr'ii- Jr Jr'ti' \PR - Pr'fi' I 

Here Sr, and Sr'^' are the cross sectional areas of the junctions in the plaquette 
R and the branches r',fj,', respectively ( the integrals are along the links between 
superconducting islands). This expression for F makes explicit its dependence on 
the geometrical characteristics of the array and the junctions. Using Eq (Eq), the 



representation of L(R, R') based on Eq (29), gives 



TIU J}l\ ^0 ^ / / dlR-dlR,dSRdSR> 

L{R,R)^ — d> f ^- . (30) 

An SrSr/ Jr Jr, \pr - PR' I 

This expression depends on the particular shape of the junctions and geometry of 
the JJA. Here we note that the general qualitative properties of the response of the 
array to external probes will not depend on the detailed form of the full inductance 
matrix. The important properties of L{R, R') are that L{R^ R') is positive, its 
non-diagonal elements are negative and decrease rapidly at long distances, plus the 
condition X]'r'=-oo ^{R^ ^) = must be fulfilled (because of the continuity of the 
flux lines). 

Here we will assume that the JJA has square lattice geometry and that each 
bond is made of cylindrical wires: Then one can obtain L{R, R'), in principle, from 
a direct integration of Eq (|30|). This is not exactly the geometry of the arrays 
considered experimentally, and we will discuss how this approximation may affect 
the final results. 

First we note that for large distances we can approximate the lattice problem 
by its continuum limit leading to 

4 
L(i?,i?')|fl-i?'|»a ~ -T-|R_ pn3 ' (31) 



which corresponds to the field of a three-dimensional magnetic dipole produced by 
a current loop. Note that the long range behavior is always given by Eq. (^) (since 
it is independent of the particular shape of the junctions). 

The short distance properties of the L{R, R') matrix depend more explicitly 
on the specific geometry of the junctions. Here we take the results of an explicit 
asymptotic evaluation of L(i?, it!')_qhi£iincd by direct integration of Eq. ( pO|) for a 
square network of cylindrical wiresE^'EZl. For example, the local behavior of L{R, R') 
is found to be given by 



2lT 



'1- KIWI) +8^-14 



L{R,R) = Lo 
L{R,R±p) = -M ^ -iLo + ^0.141875 (32) 

L{R,R±e.,±ey) = -Mn = -^0.4..., 
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where a is the lattice constant and r is the radius of the wires. For the JJA we take 
2r as the typical width of the junctions. From now on we normalize the inductance 
matrix by n^a, or A{R,R') — L{R,Fi!)/ fi^a, and Aq — Lo/^Qa, and M = M/fiaa. 
We use a/r = 10, which is a typical value for arrays made with SNS junctionsEj, 
for which, Ao = 1.13..., M = 0.14..., Mn = 0.064... . 

We have considered three different models for the inductance matrix: 
Model A is the simplest approximation that incorporates screening effects in- 
cluding only the diagonal or self-inductance contribution to L{R, R'), i.e. 

L(i?,i?') = AoMoa<5fl,fl'. (33) 

This approximation leads to null edge magnetic fie\c\&^E{R) = for all R. Model A 
has been used in thn^past by Nakajima and Sawadac3 to study vortex motion, and 
by Majhofer et al e2I to model irreversible properties of ceramic superconductors. 
This model is good whem trying to describe the properties of bulk samples, and 
three-dimensional array^. 

Model B improves model A in that it includes the nearest neighbor mutual 
inductance contributions: 

L(R, R') = Ao^ioaSR^R' - M^ioa5R,R'±^- (34) 

In this case the edge-fields as defined above are explicitly included. The corre- 
sponding magnetic fields at the boundaries are given by £'(left boundary) — MfiQa, 
i?(right boundary) = —Mhqq, with E{R) — otherwise. This model was intro- 
duced by us tostudy the efi^ect of screening currents on the nucleation of giant 
Shapiro Steps c3 and magnetic properties of JJA E3. 

Model C includes the full-range inductance matrix in the calculations. We 
model L{R, R') in this case by taking into account the local geometry of the junc- 
tions in the diagonal and mutual inductance contributions, as specified in Eq(p3), 
using the filamentary wire approximation for the remaining terms. We model E{R) 
in this case by an evaluation of Eq. (Eq) in the filamentary approximation. The 
inclusion of the full-raiige inductance matrix in the study of JJA was first consid- 
ered by Phillips et al cB, who develonfid an efficient algorithm to study static vortex 
properties in JJA. Also Reinel et a/.EH recently implemented a full inductance ma- 
trix approach to calculate static properties, somewhat closer in spirit to the method 
we iirmlemented to study the dynamics response of the JJA. More recently Phillips 
et alSB have extended their static approach to study the dynamic response of JJA. 
We will compare our results for models A through C to theirs where m)propriate. 
We have also studied the dynamical magnetic properties of Model CcJ. Here we 
concentrate on the study of the coherent vortex oscillating states of model C, which 
lead to screening-induced fractional giant Shapiro steps. 

Note that the essential difference between models A,B and C is that the former 
assumes that the magnetic field lines are constrained to lie on the plane while the 
latter takes fully into account the three dimensionality of the physical problem. Of 
course, the complexity of the algorithms needed to solve the dynamical equations 
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grows with the range of the inductance matrix. We do not discuss the specific 
details of the implementation of these algorithms, which can be found in Ref. tJ, 
but will instead concentrate on discussing the corresponding results for models A 
through C. 

4.3. Linearized equation results 

Before we discuss the numerical results obtained for models A through C we note 
that there is useful physical information that can be extracted from the linearized 
approximation to the JJA equations. Including screening current effects leads to 
a finite London penetration depth A^ . The linear approximation of the Josephson 
term in the stationary limit of Eq.( ^3) leads to the London equation for $(i?) in a 
square lattice EjO. In this case one can define an effective penetration depth that 
depends on the specific model for L{R,R'). In the local approximation of models 
A and B the penetration depth is 



with A an effective inductance, A = Aq in model A and A = Aq — 4A^ in model B. 
Since the lattice constant a plays the role of a coherence length in the JJA, we can 
define a Kj-parameter 



Kj — — — \ =- — . — =-. (36) 

a y 27r/o/ioaA y vqK 

For finite range inductance matrices, one can calculate a demagnetization factor d 
in the array. Here d is defined by the difference between the external and internal 
magnetic fields Aif — ~dM with d = ~ ^(b m J^w^R-^i^i^')' ^^^ ^^^ total 
magnetization M is proportional to J2r J{R)- In an external field the array behaves 
as if having an effective self-inductance A = (1 — d)Ao. In the diagonal case d —SL 
while in model B, d — 4A^/Ao. In the thermodynamic limit d ^ 1. One can showcJ 
that, within the linearized local approximations, for models A and B the effective 
vortex-vortex interaction energy (neglecting the lattice pinning potential) goes like 

Z^(p)=7r^j[ln— +ln2-7] for p < Aj, (37) 

P 

with 7 Euler's constant and 



e 



U{p) = rrEj-j=L== for p » Aj (38) 

y^2Trp/Xj 

with p = \R — R'\ the distance between two vortices, and the Josephson energy 

As was found originally by Pearl CJ for thin films, when considering model C, 
the long range screening currents affect the effective penetration depth and the long 
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getsEl£ 



ranffe vortex- vortex interactions. From linearizing Eq.( E7^, in the static limit, one 

27r/o^o 



Ap = t4^ (39) 



with the corresponding K-parameter, 

,^^^ = ^^^^^^_ (40) 

a ZTTlQlJQa Vg 

In this case the effective vortex- vortex interaction energy isEj, 

Z^(p)=7r£:j[ln^+ln2-7] for p < Ap, (41) 

while 

U{p) = nEj^ for p > Ap. (42) 

Here we see that U{p) has a slower decrease with p than the one given in Eq.( pq ). 
This p dependence was found in the static numericaljialculations of Phillips et al. 
E3, and it is equivalent to the one found in thin films E3. 

We see then that the main difference between models A and B, in the way we 
have defined them, is that model A does not include edge magnetic fields and model 
B does, but both have the same long range vortex-vortex interactions. In contrast, 
although models B and C both include edge magnetic fields, their corresponding 
long range vortex- vortex interactions are different. 

As we shall discuss further in this review, by considering models A through 
C we can study the effects of having local screening, long-range screening, and 
antisymmetric edge magnetic fields independently. Of course in real experimental 
samples, in principle, it is the full inductance matrix plus the antisymmetric ec 
magnetic fields that we should take into account as emphasized by Phillips et al. 

In this review we concentrate on the dynamic properties of inductive JJA. It is 
of interest, however, to mention here the results we have also obtained for the static 
magnetic and transport properties of inductive JJA as a function of k (either kj 
or Kp) for models A through C E3'E3. We can separate the results between extreme 
Type I (k <^ 1) and Type II (k ^ 1) regimes. There is a crossover k = Kx{'^ 1) 
that separates these two extreme regimes. The k > k^ regime shows a Meissner-like 
state for / < /d, with the characteristic field /d < 1, and it has a vortex lattice 
for f > fci- The vortex lattices resemble the ones that appear in the ground states 
of the extreme Type II regime (k — > oo). For k < k^ there is Type I-like behavior. 
For / < fcQ, with fco > 1, the ground state is a Meissner state. For larger fields 
/ > /co, the field penetrates the array from the boundaries of the array in the 
form of "vortex c 
form the collars t 



•s" with constant vorticity The vortices attract each other to 
. We stress that the qualitative physical aspects of the two 
extreme regimes are independent of the specific inductance model considered. The 
difference appears in the form of the long-range vortex-vortex interaction which is 
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generally only relevant for low vortex concentrations. The quantitative differences 
arise mainly in the specific values of the critical fields /d, fco, and the crossover 
screening parameter k^; E3. In the study of the subharmonic giant Shapiro steps the 
two physically extreme regimes lead to fractional giant steps in the IV characteristics 
with completely different underlying vortex dynamics. 

5. GIANT SUBHARMONIC SHAPIRO STEPS WITH SCREENING 

5.1. Coherent vortex oscillating states with local inductance matrices 

In this subsection we present the results obtained from simulations of a periodic 
JJA described by Eq. (^8|), with inductance matrices given by models A or B. Here 
we concentrate in the zero field case (/ = 0) for the same parameters used in Sec. 

3.2. In this case the dynamical equations were integrated using a second order 
Runge-Kutta method with fixed step At = 0.02 — O.lr, with r the smallest of the 
characteristic times Tg(= l/ve) or t$(= 1/v^). After discarding the results from the 
first 100 periods of the ac current, time averages were carried out for time intervals 
of, typically, up to 500 periods. 

The study of model A with / = leads to results entirely equivalent to the one 
junction case, i.e. the IV characteristics present only integer giant Shapiro steps. 
Significant changes occur, however, when rational frustrations are applied to the 
array (see Section 5.3). 

The situation for model B with / = is completely different. We studied the 
response of the inductive JJA as a function of kj, fixing AI/Aq — 0.1 (i.e. d = 0.4), 
and we found that for all the values of kj considered there are half-integer GSS. 
This appears to happen as soon as A^ ^ (i.e. whenever E{R) ^ 0). Therefore, 
the presence of the edge magnetic field term is necessary to induce subharmonic 
steps in the IVs. From an analysis of a visual animation of the vortex dynamics, 
we found that the edge magnetic fields are the basic source of VAPs that are then 
responsible for the existence of 1/2-steps in the IV characteristics. 

In Fig. 8 we show the half-integer steps in the IV curves for model B in the two 
different regimes of kj. Fig 8(a) has kj = 2.27, which corresponds to a moderate 
Type II regime. The results in Fig 8(a) are similar to the ACVS results, obtained 
when K = oo, in that they are hysteretic and the 1/2-step width is about the same 
size in both cases, for the same parameter values. 

Fig. 8(b) shows our I vs V results when kj = 0.44, i.e. in Type I regime. Again, 
there is a 1/2-step but with a hysteresis loop significantly reduced as compared to 
the one for kj — 2.21 . Furthermore, there is clear evidence of a 2/3-step in the IV, 
but with larger step width than in the k = co case and for a smaller lattice size 
(40 X 40). 

In Fig. 9 we show the change in the 1/2-step width (/ = 0) as a function of 
1/kj. The Kj = oo point, denoted as ACVS, was generated by including the edge 
magnetic fields as boundary conditions in the RSJ JJA equations (See Sec 3.5). For 
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finite kj we can distinguish three quahtatively different regimes in Fig. 9. (i) the 
Type II regime, where the 1/2-step width remains essentially constant and equal to 
the K = oo result, (ii) the Type I regime where the step width first increases and 
then (iii) decreases. The decrease occurs when the characteristic frequency v^ < v. 
This means that the current distributions are not able to follow the rapid oscillations 
due to the external current drive, and thus the coherent vortex oscillating state is 
less stable. 

We now move to discuss the microscopic coherent vortex patterns responsible 
for the 1/2-steps in model B as a function of kj. Typical examples of the family 
of oscillating coherent vortex states found in our studies of model B are shown 
in Fig. 10. In Fig. 10(a) we show results for kj — 2.27 (/ = 0). The oscillating 
vortex configurations consist of columns of mostly isolated unit charge positive and 
negative vortices with period 2/i/. As kj decreases from kj = oo to a finite value, 
1 < Kj << oo, the ACVS angle for the symmetry axis changes until it becomes 
collinear with the direction of the external current (compare with Fig. 6(a)). When 
Kj < 1, there are two distinct 1/2-step vortex oscillatory patterns one for v < v^ 
(ii) and the other for v > Vi^, (iii). In case (ii) the vortex columns are generated at 
the sample edges and then move towards the center of the array where they collide 
creating and annihilating individual vortices that concentrate about the center of 
the array. An instantaneous vortex pattern for kj — 0.44 is shown in Fig. 10(b). In 
the regime where i' > i/$, the induced currents can not follow the applied current. 
Then the vortex columns, nucleated at the sample edges, do not have enough time 
to collide with each other while still oscillating back and forth with period two. 
This situation is shown in Fig. 10(c) for kj = 0.29. In all cases considered, during 
the transient regime, VAPs are nucleated at the boundaries due to the presence 
of the antisymmetric edge-fields. But they do not correspond to the nucleation 
of commensurate vortex states as was suggested in the qualitative explanation of 
Ref. §. 

The difference in the coherent vortex patterns from the ACVS of Fig. 6(a) to 
the vortex stripes of Figs. 10(c) is due to the change in the long-range interaction 
between vortices, purely logarithmic in the non inductive ACVS to the exponential 
decay in model B given in Eq( pq) , with the Kj-dependent characteristic decay 
length from Eq( ^9|) . Of significant importance is the sign of the vortex interactions 
going from being repulsive foiJtj = oo to attractive for kj < \ (when compared to 
the lattice pinning potential 123). 

5.2. Coherent vortex oscillating states with full inductance matrix 

In this section we discuss our results from including the full-range inductance 
matrix in the analysis, within our model C approximation. This is a problem that 
has also recently been studied by Phillips et al. Ej, and we will make comparisons 
between our results and theirs where appropriate. We are also interested in com- 
paring the results discussed in this section with those of the previous section. In 
particular we want to find out how important is the form of the vortex interaction 
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at large distances, which is where the differences between model B and C arise, in 
the formation of coherent oscillating vortex patterns in zero magnetic field. 

In Fig. 11 we show our results including the full inductance matrix of model 
C, with Kp ~ 2 and for a 24 x 24 lattice. The driving ac current and frequency 
are the same as in the previous cases. The numerical integration of the dynamical 
equations is done following the same approach with time step of At = O.lrg. In 
this case we discarded the first 20 periods of the driving current and carried out 
time averages over 80 periods. In Fig. 11 we see a clearly developed hysteretic half- 
integer Shapiro step. If we compare this IV with the one shown in Fig. 8(a), for 
Kj = 2.27 and only nearest neighbours inductance, we see that the 1/2-step width 
is larger in this case. As mentioned before, the calculations with the full inductance 
matrix are more demanding than those of model B and thus we were not able to 
do a comprehensive study of the IV steps as a function of Kp, as was done in Fig. 9 
for model B. 

The important aspects of Fig. 11 are that: there is clearly a 1/2-step, there is 
hysteresis, although its size is intermediate between the kj 3> 1 and kj < 1 cases, 
and the non hysteretic part of the 1/2-step is larger than in model B. These results 
indicate that the presence of antisymmetric edge-fields {E{R) ^ 0), either modelled 
by B or C, is responsible for nucleating the coherent vortex states that lead to 1/2- 
steps in the IV's. To further check the relevance of the edge fields, we set E[R) — 
in Eq. (28) and calculated the corresponding IV. The results for model C are shown 
in Fig. 11 with a dotted line, using the same Kp and current parameters as before. 
We clearly see no evidence for a 1/2-stcp. Instead, the IV curve corresponds to the 
single junction solution times Ny. 

There are differences between models B and C in the specific vortex patterns 
formed, as we can see in Fig. 12. This is not only because the long-range interactions 
between vortices are different, but also because for model C the antisymmetric edge- 
fields are non-zero throughout the sample (in which case the strict meaning of "edge- 
fields" loses its value.) Fig. 12(a) shows the rows of isolated vortex and antivortex 
patterns for Kp — 20, as in model B. This coherent vortex state oscillates with 
period 2/j/, with vortices and antivorticcs interchanging their respective positions 
after each period T = 1/v. Note that the vortices have a wavy pattern with a 
tendency to form an angle remarkably resembling the one existing in the ACVS 
at K = oo. In Fig. 12(b) we show the vortex pattern for Kp = 2. We see that 
there is one row of vortices and one of antivorticcs, which oscillate with period 
2/v. After one period T, the rows move back and forth towards the center of the 
array, but they do not interchange positions. In Fig. 12(c) we show the vortex 
state for Kp = 0.1. In this case the large screening makes the E{R) contribution 
dominant (the magnitude of E grows as 1/Kp). In this figure we see that there are 
curved edge-field induced rows of vortices (from the left) and antivorticcs (from the 
right). Only in the center of the array can we distinguish some isolated vortices 
and antivorticcs. This coherent vortex state also oscillates with period 2/^. After 
one period T = l/iy, the edge field induced vortex rows oscillate back and forth 
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(without changing sign), whereas the isolated vortices and antivortices close to the 
center interchange their respective positions. This behavior is similar to the one 
seen in Fig. 10(b), but here the edge-field induced vortex rows occupy most of the 
sample. However, since the edge magnetic fields are size dependent, we expect that 
in real arrays (which are 1000 x 1000), the central region with the oscillating vortices 
will be bigger. 

5.3. Fractional giant Shapiro steps with screening 

We now proceed to discuss the effects of screening in the stability of the field 
induced fractional giant Shapiro steps for / = p/q Cj. In Sec. 2.3 we discussed the 
K = cxD case where the driven JJA shows giant fractional voltage steps given in 
Eq( |l^. The existence of the FGSS is due to the collective osdJlatiDn between the 
different ground state vortex configurations with frequency i^/qlS'LHala. The question 
that naturally arises is what will happen when t h ere is large screening, i.e. k < 1, 
for the ground state now is a Meissner state E3'E^. 

In this case, as we show below, we found that it was enough to use the diagonal 
inductance, model A, to simulate the JJA dynamics. The fields considered were 
/ = 1/3 and 1/2. We note that the external field induces high vortex concentrations, 
that in turn make the effect of antisymmetric edge-fields and differences in the long- 
range vortex interactions less relevant. In Figs. 13(a) and 13(b) we show the IV 
results for model A, for a lattice with 40 x 40 sites and / — 1/3. Fig. 13(a) shows the 
results for kj = 2.27, while Fig. 13(b) has kj — 0.44. It is clear from these results 
that in both limits the IV's show giant Shapiro steps at 1/3 and 1/2 fractions. We 
note that in the kj > 1 regime there is a small 1/2-step, asjvas found in the k = oo 
case using free end boundary conditions (Lee and Stroud Q). For kj < 1 the 1/2- 
step is much bigger. There are also higher harmonics seen in Fig. 13(b) that we 
have not highlighted in the figure for clarity. Fig. 13(c) shows the results for model 
B for the same parameters as in Fig. 13(b) and with d — 0.4. A comparison between 
Figs. 13(b) and (c) shows that qualitatively there are no significant differences. 

Fig. 14 shows the widths of the FGSS for / = 1/3 for steps at 1/3, 2/3 and 1/2, 
plotted as a function of 1/kj. For comparison, in Fig. 14 we also show the results 
for the 1/2-step width with / — 1/2 as a function of kj. Again, we distinguish 
three qualitatively different Kj-regimes in this figure, (i) The Type II regime where 
the step width is constant and equal to the kj = oo result; the Type I regime where 
the step width (ii) increases initially and then (iii) decreases. The decrease also 
occurs when the characteristic frequency v^ < v. As mentioned before the FGSS 
in the k = oo limit were explained in terms^^ coherent oscillation of the ground 
state vortex lattices formed when J"^^* = ErLlu'EI. When screening is included the 
vortex lattice configurations get modified. In Fig. 15 we show the instantaneous 
vortex configurations at the 1/2-step in the / = 1/2 case, both in the Type II and 
Type I cases. In Fig. 15(a) we took kj — 2.27, and the vortex state shown in 
the figure oscillates with period 2/i/. We notice that there are remanents of the 
ground states present in the k = oo limit but now forming "clusters" separated by 
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"Bloch" or "soliton" surface walls due to the screening currents. The size of the 
"cluster" regions in this particular case is about £ = 7a while the Bloch wall is 2a. 
As Kj —^ Kx, with Kj; < 1 the number of Bloch walls increases, and for kj < Kx we 
notice the merging of vortices into the stripes shown in Fig. 15(b). For kj = 0.44. 
This structure is reminiscent of the intermediate state in a bulk superconductor, 
although the structure discussed here is generated dynamically. The 1/2-step in 
the IV-characteristics corresponds to a coherent oscillation of this stripped vortex 
distribution that oscillates between +1 and zero vorticities with period 2 /v. In the 
Type I regime with J'^^* = . the equilibrium magnetic field distributions correspond 
to the Meissner state E3cJ. Starting from this state and turning on the external 
I'^^^it) current, the induced Lorentz force "pulls" the flux into the sample and, 
after a transient, the vortex distributions settles into the final oscillating stripped 
configuration. This is shown in Fig. 16 for / = 1/3, within the 1/3-step. 

In both Type I and Type II limits, the oscillating vortex states are non-equilibrium 
stationary states. In the Type II regime the oscillating vortex lattice, even when 
similar to the ground state of the k j = oo case, is different because of the screening 
induced defects. The vortex state without external driving, /'^^* = 0, forms a lattice 
but with a density of vortices smaller than_t]p£_external field density / = p/q, be- 
cause there is now a finite penetration depth E3'E^. On the other hand, the oscillating 
vortex state at the fractional steps, has a vortex density exactly equal to f ~ p/q, 
so that it can oscillate with frequency I'/q. Only in the k = oo limit both the oscil- 
lating vortex state and the ground state coincide. In the Type I regime, as it was 
mentioned before, the difference is more explicit. Whereas the ground state is the 
Meissner state, without vortices, the oscillating vortex state is the "intermediate" 
state shown in Fig. 15(b) and Fig. 16(d). 

Although we did not carry out calculations for / = p/q with model C we do 
not expect significant qualitative changes as compared to the picture presented 
above, since for high vortex concentrations we have not found significant differences 
between models A and C in the static case l3. 

5.4. Magnetic field dependence of half-integer Shapiro steps 

In the experiments by H. C. Lee et al. Ej it was found that the 1/2-steps exist 
not only for zero field, but also for any finite value of /. As we discussed above, we 
have found that as long as we include the antisymmetric edge-fields, represented by 
a E{R) 7^ 0, there are 1/2-steps in the IV's for all the values of / considered. In 
Fig. 17 we show the step width of the 1/2-integer step as a function of / both for a 
Type II [kj — 2.27) and Type I cases {kj = 0.44). In this case we simulated model 
B with d = 0.4. In the Type II case the curve shows more structure than in the 
Type I regime. Also, in the Type II regime the 1/2-step width close to / = tends 
to go to zero for small fields. Moreover, when the ACVS is induced by disorder 
{k = oo) we found that the 1/2-step disappears for fields smaller than / = 0.05. 
What is qualitatively important is that in both regini£s the curves have maxima at 
/ = and / = 0.5, as observed in the experiments Ej. The relative height of the 
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two maxima depends on the frequency and amplitude of the ac current. 

Ahhough all the results presented above are T — Q results, we have verified that 
they are stable at sufficiently low temperatures. 

6. CONCLUDING REMARKS 

The central theme of this review was the search for a microscopic understand- 
ing of the different types of steps that can be found in the IV characteristics of 
driven proximity effect Josephson junction arrays (JJA). The motivation for this 
work started with the experimental findings of giant integer (/ = 0) and fractional 
(/ = P/l) Shapiro steps BB LGSS), that were theoretically explained in terms of an 
extreme Type II RS J model Q'Bcl. The dynamical equations of motion describing the 
JJA are equivalent to a large set of overdamped nonlinear coupled equations driven 
by an external time-dependent current. However, in the / = case the array be- 
haves as a set of single junctions connected in series, and thus the GSS are simply a 
manifestation of the Shapiro steps of each junction. Note that if we were discussing 
the underdamped limit of a single junction extra subharmonic states would appear. 
When / — p/q the steps seen in the IV's are due to a coherent oscillation of the 
ground state vortex lattice B'LI. 

Experimentally, there are fractional giant Shapiro steps in zero field that do 
not conform to the picture given aboveOO'EJO. In order to shed light into this 
problem we decided to study the effects of disorder in a driven JJA modeled by 
the RSJ model in zero field and k = oo ESiIj'Lj. We discovered that as soon as 
we break translational invariance the corresponding IV's show half-integer GSS. 
When looking in more detail at the corresponding microscopic vortex dynamics, we 
found a stationary coherent oscillatory vortex pattern with completely unexpected 
properties. We termed this novel state axisymmetric coherent vortex state (ACVS), 
which characterizes its geometric nature. We have carried out a comprehensive 
analysis of the stability properties of the ACVS against variations of several physical 
parameters. We have concluded that the ACVS is in fact a very stable global 
attractor of the nonlinear dynamics, that can be generated in several different ways. 
The only essential element needed to nucleate the ACVS, however, is to have a 
mechanism that produces a vortex antivortex pair (VAP) in the initial conditions. 

Unfortunately, we still do not have a complete theoretical understanding of the 
ACVS, in particular, the "radiation" and annihilation mechanisms that are essential 
for its formation and stability. The VAP radiation phenomenon, due to the presence 
of the ac current, does not appear to have been studied before. Furthermore, we 
would like to understand why the ACVS forms in the very specific patterns shown 
in the figures. Note that the ACVS is a true nonequilibrium state that disappears 
as soon as we turn off both the dc and/or ac currents. This property seems to imply 
that the ACVS can not be explained from a minimum energy principle. This is a 
fundamental difference between the ACVS and the fractional GSS observed with 

/ = p/qm 

Most of the experiments that show GSS were made in proximity effect arrays, 
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with large screening currents. We began the study of self- field effects in JJA by 
considering local screening models for the inductance niatrixcJ. There we found 
that there are new fractional GSS triggered by the current induced antisymmetric 
magnetic edge-fields. We found that there is a robust 1/2-step in the fV's even 
when / = 0, and similar to what was seen in the experiments. 

Recently, Phillips et al. Ej and Reinel et al. l3 have extended the static analysis 
including the full inductance matrix, and then Phillips et al. c3 extended our 
dynamic study including the full inductance matrix. The inductance matrix model 
used by Phillips et al. takes fully into account the specific geometry of typical 
SIS arrays. The local geometry of these arravs ia^ quantitatively different from the 
SNS junctions used in the study of GSS HEjiij^EII. However, their results must be 
semi-quantitatively correct, in particular at long distances. 

In this review we have presented our own full inductance matrix results within 
the filamentary approximation. We have discussed the main differences between 
local screening models (A and B) and the full inductance model (C). Models A 
and B differ in that model A does not include "edge-fields" {E{R) = 0) while 
model B does {E{R) ^ 0). Nonetheless, as shown in Section 4.3, both models 
have, within the linearized approximation of the equations of motion, the same 
long distance vortex-vortex (V-V) interaction potential. On the other hand, models 
B and C have different long distance V-V interaction potential. In model B the 
V-V interaction decays exponentially with distance, while in model C it decays 
algebraically. Nevertheless, both models include the "edge-fields" (quantitatively 
the range of E{K) depends on the range of L{R,R').) 

From our comparisons of the three inductance matrix models, the following 
picture emerges: 

(i) The existence of the subharmonic GSS at / = is due to the presence 
of current induced antisymmetric magnetic fields, which break the translational 
symmetry. In fact, when we artificially imposed E{R) = in model C, we did 
not obtain a subharmonic response, despite of the fact that the long range V-V 
interactions were properly taken into account. Instead, we get the single junction 
type IV (see Fig. 11.) 

(ii) The specific structure of the underlying coherent oscillating vortex states and 
their corresponding subharmonic step widths, depend directly on the long range V- 
V interactions (as also found by Phillips et al..) We found that as the screening 
parameter k decreases, there is a whole family of coherent oscillating vortex states 
where the vortex rows tend to line up with the external current, for a given in- 
ductance model. For k larger than the crossover Kx{^ 1), the V-V interaction is 
repulsive while for k < k^ it is attractive. This explains why the vortex patterns in 
the former case contain isolated vortices while in the latter the vortices form stripes 
of constant vorticity. Also, the geometric vortex patterns depend on the model used 
to describe the inductance matrix for a given k, which is connected to the difference 
in the nature of the V-V interaction. 

In all the cases studied (either with f — or f = p/q) the coherent vortex 
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states arc a fundamentally non-equilibrium consequence of the dynamics, with the 
only exception being the particular case of k = c» and / = p/q- We believe that 
these non-equilibrium coherent vortex states, due to the breaking of translational 
invariance, are the underlying mechanism behind all the fractional giant Shapiro 
steps observed in two-dimensional SNS JJA. At / = p/ q this translational symmetry 
is broken by the applied magnetic field, that generates higher order periodicity thus 
leading to fractional GSS. At / = the translational symmetry can be broken 
either by defects or by the current induced antisymmetric magnetic fields. They 
tend to nucleate vortex-antivortex pairs that produce ACVS-like vortex patterns 
responsible for the subharuaonic GSS. Fractional steps can also be generated in 
lower dimensional systems by breaking translational invariance (e.g. ladder arrays 
and single plaquettes EaoO'Cj) . 

The subharmonic GSS studied here correspond to fully overdampcd JJA. It is 
known that underdanmed single Josephson junctions can have subharmonic Shapiro 
steps and even chaos e2I. In fact^^ome chaotic phenomena has been recently studied 
in ac driven underdamped JJA e3. It would be interesting also to know the effect of 
a finite capacitance (but not large enough to induce chaos) in these coherent vortex 
states, together with a finite vortex mass and the effects of the non-linear vortex 
viscosity c3. 

We^lso want to mention other possiblyjelated systems, like charge density 
waves e3 and the Frenkel-Kontorova model Eil, where subharmonic Shapiro steps 
have been found as a consequence of the collective nonlinear dynamics of these 
systems. 

How can the presence of a specific coherent oscillatory vortex state be verified 
experimentally? There are several possibilities but we will only mention the ones 
that are closest to the JJA system. Recently, many different techniques have been 
developed for doing spatially resolved measurements in JJA oEjOH (save for 
Tonomura's al. approacho.) Most of these techniques, however, are dedicated 
to the study of static flux configurations E§E3l£3. Lachenmann et al cB have been 
able to study dynamic states in JJA by the using the low temperature scanning 
electron microscopy (LTSEM) technique. These are basically measurements of the 
local distribution of dissipation produced by the local heating of the LTSEM tip, 
that in turn gives the average junction voltage in the arrays. Since a moving vortex 
generates a voltage as it crosses a junction, this type of measurement indicates 
the location of the vortices during their oscillations (note that since this process 
is independent of the sign of the vorticity, it will not average to zero after many 
ACVS ac current periods.) In fact, in Fig. 18, we show a simulation of the average 
voltage along the current direction for the ACVS at a finite temperature. There 
we see where the vortices were located most of the time regardless of their sign. 
In particular, the characteristic angle of the ACVS is clear. Therefore, an LTSEM 
based measurement may be able to provide evidence for the formation of an ACVS. 

In conclusion, we have shown that giant subharmonic steps in the IV charac- 
teristics can be generated under many different circumstances at zero and non-zero 
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external magnetic field, but for which the underlying microscopic vortex dynam- 
ics can be fundamentally different. In all the cases considered, at least for two 
dimensional SNS arrays, the steps appear to always entail a microscopic coherent 
oscillating vortex state. 

We do not as yet have a more theoretical understanding of all the nonlinear 
mechanisms responsible for the formation of all these coherent vortex states and 
further studies are needed to shed light on these questions. 
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7. Figure Captions 



Figure 1: Schematic representation of a Josephson junction square array driven 
by an external current /. Tlie sites denote superconducting islands and the crosses 
the junctions themselves. 



Figure 2: Vortex configurations, for a 40x40 array with one defect, produced 
by a dc current with / = fO, 5x ~ 0.1, Sy = 0. The white and black squares 
denote —1 and +1 phase vorticity per plaquette (see Eq. (11)), respectively. No 
squares here means zero vorticity. For currents (a) I^c = O.S/q < Ic ^ 0.84/o; (b) 
Idc = 0.87/0 > Ic- 



Figure 3: IV characteristics for a 40 x 40 square array with lac = Iqt v = O.lz^o 
and pbc. (a) Ordered array; (b) array with one defect in the center of the lattice 
and 5x — 0.1, 5y — and / = 2 (for clarity the IV curve is shifted by one unit); (c) 
blowup of (b) for an 80 x 80 array showing hysteretic behavior and a clear 2/3-step. 



Figure 4: Step width for of the 1/2-stcp for arrays with one defect given by the 
same parameters as in Fig. 3(b). (a) As a function of the frequency, {lac ~ Iq)', and 
(b) as a function of the ac current, (v = 0.1 vq). 



Figure 5: Number of vortices as a function of time for an array with one defect 
and I(ic = 0.586/o (within the 1/2-integer step). The lattice size, frequency and ac 
current are the same as in Fig. 3 (b). (a) Absolute number of vortices Na{t); (b) 
total number of vortices Nxit). See text for definitions of td, tAcvs and toff- Here 
time is measured in units of l/i^o- 



Figure 6: Axisymmetric coherent vortex states (ACVS). (a-b) Stationary oscil- 
lating vortex configuration for an ACVS in the 1/2-step, for the same parameter 
values as in Fig. 3(b) with I^c = 0.586/o. (c) ACVS configuration in the 2/3-stcp 
with Idc = 0.622/o in an 80 x 80 lattice. We use the same conventions as in Fig. 2. 



Figure 7: Vortex distributions showing an ACVS at a 1/2-integer step generated 
by antisymmetric edge-fields in a 64 x 64 lattice with 7=1, lac = Iot '^ = O.lz^o 
and Idc — 0.57/0. The notation is the same as in Fig. 2. 
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Figure 8: IV-characteristics for arrays with only nearest-neighbours inductance 
(model B): Here / = 0, h.c. = h, v ^ O.lvg, and M = O.IAq {d = 0.4), for 40 x 40 
lattice size, (a) Type II array with kj = 2.27; (b) Type I array with kj = 0.44. 



Figure 9: IV-step widths AI/Iq vs 1/kj, for half-integer steps in zero field. 1^ , 
Io,iy = O.We,M = O.lAo, and lattice size 40 x 40. 



Figure 10: Instantaneous vortex configurations for 1/2-steps with / = (model 
B). Same A4 and ac driving currents as in Fig. 9, and Idc = O.6I/0. (a) kj = 2.27; 
(b) Kj — 0.44; and (c) kj ~ 0.29. We use the same conventions as in Fig. 2 to 
denote vorticity. 



Figure 11: IV-characteristics of an array with full inductance matrix (model C): 
Kp = 2, lattice size 24 x 24, / = 0, la.c. = Iq, and v = O.li^g. The dotted line 
corresponds to a simulation with the same model C inductance, same Kp and applied 
currents, but without edge fields (i.e. imposing E{R) = in Eq.(24)). 



Figure 12: Instantaneous vortex configurations for 1/2-steps at / ~ 0, lac ~ lo, 
V — O.li'e, and Idc = 0.59/o, for model C. (a) Kp = 20; (b) Kp = 2, (c) Kp = 0.1. 
The same conventions to indicate vorticity as in Fig. 2 is used. 



Figure 13: IV-characteristics for model A with / = 1/3, la.c. = lo, and jy = O.li^g, 
for lattice size 40 x 40. (a) Type II array with kj = 2.27; (b) Type I array with 
Kj — 0.44; (c) The same as (b) including edge-fields (model B) with Al/Ag = 0.1 
(d = 0.4). The giant Shapiro steps at n/3 and n/2 are marked. Curves start about 
the same place and are displaced from each other for clarity. 



Figure 14: IV-step widths, AI/Ic vs 1/kj, for la.c. = lo,'^ = O.li'e, diagonal 
inductance approximation (model A), and lattice size 40 x 40. Upper curve corre- 
sponds to 1/2-step for / = 1/2, while the three lower curves to 2/3-step, 1/3-step 
and 1/2-step with / = 1/3. 



Figure 15: Instantaneous vortex configurations for 1/2-steps with / = 1/2, in the 
diagonal inductance approximation (model A). lac — li), v — Q.lvg. (a) kj ^ 2.27 
(Type II), with Idc = 0.4/o. (b) kj = 0.44 (Type I), with Ijc = O.6/0. White and 



34 Daniel Dominguez and Jorge V. Jose 

black squares indicate and +1 phase vorticity per plaquette, respectively. 



Figure 16: Transient time evolution for the vortex configurations at a 1/3-step for 
/ = 1/3, in model A. Here, Idc = 0.45/o, lac = lo, smd v = O.ljyg. The current 
is applied along the vertical direction, (a) Initial Meissner configuration (i.e. the 
ground state for I^xt = 0), (b) The first vortex row has penetrated the array, (c) 
Advancing front of vortex rows, (d) Final oscillating 1/3-period state of vortex 
rows. White and black squares indicate and +1 phase vorticity per plaquette, 
respectively. 



Figure 17: Current width of the 1/2-steps as a function of / (model B) for: (a) 
Kj — 2.27 and (b) kj = 0.44. Same Ai and driving currents as in Fig. 9. 



Figure 18: Time average of the voltage drop in each junction along the current 
direction. Contour plot for the case with Idc — 0.57/o, lac = loi ^ = O.liyg (i.e. 
in the half-integer step), for a 40 x 40 lattice with k = oo and finite temperature 
ksT = O.OlEj. 



